Part I

Load Plot and Shapefile

library( maptools )
library( sp )
library(maps)

syr <- readShapePoly( fn="shapefiles/tl_2010_36067_tract10", proj4string=CRS("+proj=longlat +datum=WGS84") )
syr <- syr[as.numeric(as.character(syr$NAME10)) < 64, ]
par( mar=c(0,0,0,0) )
plot( syr,  border="gray10" )

Add Roads

roads <- readShapeLines( fn="shapefiles/tl_2015_36_prisecroads", proj4string=CRS("+proj=longlat +datum=WGS84") )

#Extract Intersate
interstate <- roads[ roads$RTTYP == "I" , ]

Look at unclipped portion

par( mar=c(0,0,0,0) )
plot( syr,  border="gray80" )
plot( roads, col="steelblue", lwd=2, add=T)
plot(interstate, col = "red", add = T)

Clip

library(raster)
library(rgeos)

interstate_clipped <- gIntersection(syr, interstate)

Plot Clipped

plot( syr,  border="gray80" )
#plot( roads_clipped, col="steelblue", lwd=2, add=T)
plot( interstate_clipped, col="red", add=T)
map.scale( metric=F, ratio=F, relwidth = 0.05, cex=0.5 )

Houses

library(dplyr)
houses <- read.csv("https://raw.githubusercontent.com/lecy/hedonic-prices/master/Data/Housing%20Price%20In-Class%20Exercise%20(Responses).csv", stringsAsFactors = F)

lat.lon <- read.csv("lat.lon.csv", stringsAsFactors = F)
lat.lon <- lat.lon[ c("lon","lat")]

plot( syr,  border="gray80" )
#plot( roads_clipped, col="steelblue", lwd=2, add=T)
plot( interstate_clipped, col="red", add=T)
map.scale( metric=F, ratio=F, relwidth = 0.05, cex=0.5 )
  
points(lat.lon$lon, lat.lon$lat, pch = 19, cex=.3)
buff <- gBuffer(interstate_clipped, width = .005)
plot(buff, add=T)

Less find the points near dat highway

lat.lon <- SpatialPoints(lat.lon, proj4string = CRS("+proj=longlat +datum=WGS84") )

overStuff <- over(lat.lon, buff)
overStuff[is.na(overStuff)] <- 0
nearHighway <- overStuff==1

houses<- mutate(houses, nearHighway)

nearHighway
##     1     2     3     4     5     6     7     8     9    10    11    12 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
##    13    14    15    16    17    18    19    20    21    22    23    24 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
##    25    26    27    28    29    30    31    32    33    34    35    36 
## FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
##    37    38    39    40    41    42    43    44    45    46    47    48 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE FALSE 
##    49    50    51    52    53    54    55    56    57    58    59    60 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
##    61    62    63    64    65    66    67    68    69    70    71    72 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
##    73    74    75    76    77    78    79    80    81    82    83    84 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE 
##    85    86    87    88    89    90    91    92    93    94    95    96 
## FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
##    97    98    99   100   101   102   103   104   105   106   107   108 
## FALSE FALSE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE 
##   109   110   111   112   113   114   115   116   117   118   119   120 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
##   121   122   123   124   125   126   127   128   129   130   131   132 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
##   133   134   135   136   137   138   139   140   141   142   143   144 
## FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE 
##   145   146   147   148 
##  TRUE FALSE FALSE FALSE

Less plot it

plot(syr, border = "gray80")
plot(interstate_clipped, col = "red", add = T)
points(lat.lon$lon[!nearHighway], lat.lon$lat[!nearHighway], pch = 19, cex = .3)
points(lat.lon$lon[nearHighway], lat.lon$lat[nearHighway], pch = 19, cex = 1, col = "yellow")

Part II

library( rgdal )
library( maptools )
library( geojsonio )
library(dplyr)
library(RColorBrewer)
library(maps)

syr <- readOGR(dsn = "shapefiles/syr_parcels.geojson")
## OGR data source with driver: GeoJSON 
## Source: "shapefiles/syr_parcels.geojson", layer: "OGRGeoJSON"
## with 41502 features
## It has 64 fields
industrials <- syr[syr$LandUse=="Industrial", ]
notIndustrials <- syr[syr$LandUse!="Industrial", ]

plot(notIndustrials, col = "grey", border = F)
plot(industrials, add = T, col = "blue", border = F)

##Buff it up

plot(notIndustrials, col = "grey", border = F)
plot(industrials, add = T, col = "blue", border = F)

industrialBuff <- gBuffer(industrials, width = .003621, capStyle = "FLAT", quadsegs = 1)
plot(industrialBuff, add = T, col = rgb(44, 4, 248, alpha = 50, maxColorValue = 355))
map.scale( metric=F, ratio=F, relwidth = 0.15, cex=0.5 )

##Add houses

plot(notIndustrials, col = "grey", border = F)
plot(industrials, add = T, col = "blue", border = F)

plot(industrialBuff, add = T, col = rgb(44, 4, 248, alpha = 50, maxColorValue = 355))
map.scale( metric=F, ratio=F, relwidth = 0.15, cex=0.5 )
points(lat.lon$lon, lat.lon$lat, pch = 19, cex=.2)

Let’s get some overs

lat.lon <- read.csv("lat.lon.csv", stringsAsFactors = F)
lat.lon <- lat.lon[ c("lon","lat")]
lat.lon <- SpatialPoints(lat.lon, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") )
overs <- over(lat.lon, industrialBuff)
nearIndustry <- overs==1
names(nearIndustry) <- NULL
nearIndustry[is.na(nearIndustry)]<- F
houses <- cbind(houses, nearIndustry)

Let’s plot this

plot(notIndustrials, col = "grey", border = F)
plot(industrials, add = T, col = "blue", border = F)

points(lat.lon$lon[nearIndustry], lat.lon$lat[nearIndustry], pch = 19, cex=.4, col = "yellow")
points(lat.lon$lon[!nearIndustry], lat.lon$lat[!nearIndustry], pch = 19, cex=.2)

Now with schools

schools <- syr[syr$LandUse=="Schools", ]
notSchools <- syr[syr$LandUse!="Schools", ]

plot(notSchools, col = "grey", border = F)
plot(schools, add = T, col = "red", border = F)

##Less add the buffer

plot(notSchools, col = "grey", border = F)
plot(schools, add = T, col = "red", border = F)

schoolBuff <- gBuffer(schools, width = .001812, capStyle = "FLAT", quadsegs = 1)
plot(schoolBuff, add = T, col = rgb(248, 18, 18, alpha = 50, maxColorValue = 355))
map.scale( metric=F, ratio=F, relwidth = 0.15, cex=0.5 )

##Less plot some houses

plot(notSchools, col = "grey", border = F)
plot(schools, add = T, col = "red", border = F)

plot(schoolBuff, add = T, col = rgb(248, 18, 18, alpha = 50, maxColorValue = 355))
map.scale( metric=F, ratio=F, relwidth = 0.15, cex=0.5 )
points(lat.lon$lon, lat.lon$lat, pch = 19, cex=.2)

Let’s over

lat.lon <- read.csv("lat.lon.csv", stringsAsFactors = F)
lat.lon <- lat.lon[ c("lon","lat")]
lat.lon <- SpatialPoints(lat.lon, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") )
overs <- over(lat.lon, schoolBuff)
nearSchools <- overs==1
names(nearSchools) <- NULL
nearSchools[is.na(nearSchools)]<- F
houses <- cbind(houses, nearSchools)

Less plot it

plot(notSchools, col = "grey", border = F)
plot(schools, add = T, col = "red", border = F)

points(lat.lon$lon[nearSchools], lat.lon$lat[nearSchools], pch = 19, cex=.4, col = "yellow")
points(lat.lon$lon[!nearSchools], lat.lon$lat[!nearSchools], pch = 19, cex=.2)